* Settings
global SSDIMed "/disk/agedisk4/medicare.work/miller-DUA50377/proj_ssdi"
version 16
do "$SSDIMed/scripts/_auxiliary/_project_settings.do"

* Confirm that the Stata global for the project has been defined
assert !missing("$SSDIMed")

* All required Stata modules are available in the /_packages and /_auxiliary folders
cap adopath - "/home/site/etc/stata/ado.nber"
cap adopath - "/home/site/etc/stata/ado.nber/updates"
adopath ++ "$SSDIMed/scripts/_packages"
adopath ++ "$SSDIMed/scripts/_auxiliary"
capture log close 
log using $SSDIMed/scripts/logs/05_fig6c_marginal_mortality_fn.log, text replace 

cap do "$SSDIMed/scripts/_auxiliary/_project_grstyle.do" 3in 3in

************************************
***MTO for mortality 
local DeltaC=-5000
local posDeltaC=-1*`DeltaC'
***use county-month clusters or county clusters 
local clustering ""
*local clustering "_cluster_county"
use "$SSDIMed/results/estimates/marginal_function/tablecolumn2`clustering'.dta", clear 
rename (lb2 sd2 ub2 value2) (lb sd ub value)
drop sd 
drop if (regexm(var, "500") & !regexm(var, "5000")) | regexm(var, "50000")
keep if regexm(var, "x51") | regexm(var, "x51UR") | regexm(var, "x52") | regexm(var, "x52UR") | regexm(var, "00")
drop if regexm(var, "C_") | regexm(var, "diff_")
replace var=subinstr(var, "`posDeltaC'", "", 1) 
gen x=.
gen xlb=.
gen xub=.
foreach var in x51 x51UR x52 x52UR {
replace x=value if var=="`var'"
replace xlb=lb if var=="`var'"
replace xub=ub if var=="`var'"
}
replace x=200 if var=="B_200" | var=="BUR_200"
replace x=800 if var=="B_800" | var=="BUR_800"
gen B=. 
gen Blb=.
gen Bub=.
foreach var in B_x51 B_x51UR B_x52 B_x52UR B_200 BUR_200 B_800 BUR_800{
replace B=value if var=="`var'"
replace Blb=lb if var=="`var'"
replace Bub=ub if var=="`var'"
}
 gen UR=regexm(var, "UR")
replace var=subinstr(var, "B_", "", 1)
replace var=subinstr(var, "BUR_", "", 1)
keep if inlist(var, "x51", "x51UR", "x52", "x52UR", "200", "800")
 * foreach var in B Blb Bub {
 *   replace `var'=`var'/1000
 * }


***the mins and maxes below only affect the 200 & 900 that we use for the CIs 
gcollapse x xlb xub B Blb Bub , by(var UR) 

foreach var in x51 x51UR x52 x52UR{
sum x if var=="`var'" 
local `var': di %3.0f = `r(mean)'
}

* Override default graph sizes
grstyle set size 8pt: text_option
grstyle set size 9pt: axis_title
grstyle set size 8pt: tick_label
grstyle set size 2pt: tickgap

* Arrows and lines
local p 5
grstyle set color gs8: p`p'
grstyle set lpattern solid: p`p'
grstyle set linewidth vvthin: p`p'

* How many major ticks?
* How much should range extend beyond ticks (as fraction of tick step size)
local y_ticks 7
local y_margin_lo 0.4
local y_margin_hi 0.4

* y1 tick lo, step size, format
local y1_tick_lo 250
local y1_step 10
local y1_fmt "%9.0f"
local y1_sign ""
local y1_unit ""

* Construct y labels and scales from parameters above
forvalues y = 1/1 {
  local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
  local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
  local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
  di "y`y'_tick_hi:  `y`y'_tick_hi'"
  di "y`y'_range_lo: `y`y'_range_lo'"
  di "y`y'_range_hi: `y`y'_range_hi'"
  
  local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
  local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
  di `"y`y'_label:  `y`y'_label'"'
  di `"y`y'_scale:  `y`y'_scale'"'
}

* x-axis settings
local xlo 120
local xhi 830
local xscale xscale(range(`xlo' `xhi'))
local x_axis xlabel(200(200)800) `xscale' xtitle("Entrants per million", margin(t=-.5))

* style and color for y-axis title/labels
local y1_title_style it
local y1_title_color gs8
local y1_label_style it
local y1_label_color gs8
local y1_label_size 8pt

local text_style sf
local text_color black

* y-axis labels placed above grid rule (build manually)
local y1_x `xlo'
local y1_place ne
local y1_just left
local y1_width 25

local y 1
local y`y'_la
foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
  if `tick' == `y`y'_tick_hi' {
    * white textbox to go behind label
    local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(`y`y'_width') height(4) bcolor(white%70))
    *label
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
  else {
    local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
  }
}

* y-axis titles, horizontal orientation (build manually)
local y 1
local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Deaths per 10,000}"
local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))

* y-axis settings
local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'

* Plot lines
local y1 B
local xvar x
local line_y1_lfit "(lfit B x if UR==0, col(orange) lw(thick) range(200 800))"
local line_y1_scat "(scatter B x if UR==0 & inrange(x, 201, 799),msym(S) mcolor(black) msize(*.8))"
local line_y2_lfit "(lfit B x if UR==1, lpat(dash) lcol(orange) lw(thick) range(200 800))"
local line_y2_scat "(scatter B x if UR==1 & inrange(x, 201, 799),msym(Oh) mcolor(black) msize(*1) mlwidth(medthick)) "
  local ci_area_0 (rarea Blb Bub `xvar' if UR==0, sort col(gs8%30) lw(none))
  local ci_area_1 (rarea Blb Bub `xvar' if UR==1, sort col(gs8%30) lw(none))

* Label line_y1
local text_1 text(280 500 "{`text_style':Mean unemployment}"
local text_1 `text_1', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
local xref = `xvar'[4]
sum `y1' in 4 
local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + .7' `=`xref' + 4.5', pstyle(p5))

* Label line_y2
local text_2 text(280 550 "{`text_style':+1 p.p. unemployment}"
local text_2 `text_2', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
local xref = `xvar'[4]
sum `y1' in 4 
local pci_2 (pci `=r(mean)' `=`xref'' `=r(mean) + 0.07' `=`xref' + 0.45', pstyle(p5))

***11.5=240
***12=250
***12.5=260
***13=270
***13.5=280
***14=290
grstyle set graphsize 3in 3in
local legend legend(off)
local region graphregion(margin(t=5 b=-1))
local ws "placement(w) size(*1)"
local es "placement(e) size(*1)"
local labels          text(281.5 290 "{`text_style':49-year-olds}" "{`text_style':at mean}" "{`text_style':unemployment}", placement(c)) 
local pci_labels              (pci 285 362 288 375, pstyle(p5))
local labels `labels' text(295 510 "{`text_style':49-year-olds}" "{`text_style':at +1 p.p.}" "{`text_style':unemployment}", placement(c))
local pci_labels `pci_labels' (pci 291.5 422 289.5 396, pstyle(p5))
local labels `labels' text(265 510 "{`text_style':50-year-olds}" "{`text_style':at mean}" "{`text_style':unemployment}", placement(c)) 
local pci_labels `pci_labels' (pci 268 582 272.5 627, pstyle(p5))
local labels `labels' text(285 712 "{`text_style':50-year-olds}" "{`text_style':at +1 p.p.}" "{`text_style':unemployment}", placement(c))
local pci_labels `pci_labels' (pci 280.5 710 273 682, pstyle(p5))
twoway `ci_area_0' `ci_area_1' `line_y1_lfit' `line_y2_lfit' `line_y1_scat' `line_y2_scat' `pci_labels', `x_axis' `y_axis1' `legend' `region' `labels' graphregion(color(white))  bgcolor(white) 

* Export graph
cap mkdir "$SSDIMed/results/figures"
graph export "$SSDIMed/results/figures/estimatedmodel_mortality4950`clustering'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace

** EOF






